/*
文件名：p_comp.h
时间：2024-4-27
功能：
	电力信号的压缩
*/
#include "p_comp.h"

u32 com_time_getus(void)
{
	auto now=chrono::system_clock::now();
	return chrono::duration_cast<chrono::microseconds>(now.time_since_epoch()).count();
}

float f_s_in[S_N]; //浮点信号，fft输入用，后续用作残差
s16 s16_s_in[S_N]; //主峰时域波形
float tmp_fre_buf[S_N]; //调试用频域缓存
int fft_ip[S_N/2]; //用(2 + sqrt(FFT_LEN))，给足够大
float fft_w[S_N/2]; //用FFT_LEN / 2
void tc_compress(s16 *in, S_PCP_MAIN_FRE &out) //in: 输入, out: 输出缓存
{
	int i=0;
//1、做实数fft
	for(i=0;i<S_N;i++)
	{
		f_s_in[i]=in[i];
	}
	fft_ip[0]=0; //每次都重新算

	u32 t0 = com_time_getus();
	rdft(S_N, 1, f_s_in, fft_ip, fft_w); //正变换
	u32 t1 = com_time_getus();
	printf("fft: %d\n",t1-t0);
//2、取主峰
	int max_ind=1;
	float max_val=0;
	for(i=0;i<S_N/2;i++) //复数个数
	{
		int ind=i*2;
		float amp=sqrt(f_s_in[ind]*f_s_in[ind] + f_s_in[ind+1]*f_s_in[ind+1]);
		tmp_fre_buf[i]=amp;
		if(amp>max_val)
		{
			max_ind=i;
			max_val=amp;
		}
	}
	int st=max_ind-P_W2;
	int end=max_ind+P_W2; //取主峰前后P_W2个
	if(st<1) st=1;
	if(end>=S_N/2) end=S_N/2-1;

	out.main_st=st;
	int tst=st;
	for(i=0; i<ARRAY_SIZE(out.fre_buf)/2 && tst<=end; tst++,i+=2)
	{
		out.fre_buf[i*2]=f_s_in[tst*2];
		out.fre_buf[i*2+1]=f_s_in[tst*2+1];
	}
//3、反变换，得到偏差
	for(i=0;i<S_N/2;i++) //复数个数
	{
		if(i<st || i>end)
		{
			int ind=i*2;
			f_s_in[ind]=0;
			f_s_in[ind+1]=0;
		}
	}
	rdft(S_N, -1, f_s_in, fft_ip, fft_w); //反变换
	for(i=0;i<S_N;i++) //与输入的做差
	{
		f_s_in[i]=f_s_in[i]/S_N*2; //这个库做出来扩大了N/2倍
		s16_s_in[i]=f_s_in[i]; //主峰时域波形
		f_s_in[i]=in[i]-s16_s_in[i]; //与输入做差
		//f_s_in[i]/=2; //缩小幅度，有损压缩
	}
	for(i=1;i<S_N;i++) //误差的差分值
	{
		out.dd_list[i] = (f_s_in[i] - f_s_in[i-1]);
		//out.dd_list[i] = f_s_in[i]; //直接输出误差较大，压缩率下不去
	}
	out.dd_list[0]=f_s_in[0]; //第一个元素是误差初始值

//调试
	return ;
	for(i=0;i<S_N;i++) //误差的差分值
	{
		//输入数据	频域数据	主峰时域	误差差分
		int ind=i*2;
		//float main_amp=sqrt(f_s_in[ind]*f_s_in[ind] + f_s_in[ind+1]*f_s_in[ind+1]);
		printf("%d	%f	%d	%d\n",
			in[i], tmp_fre_buf[i], s16_s_in[i], (int)(out.dd_list[i]));
	}
}

void tc_decompress(S_PCP_MAIN_FRE &in, s16 *out)
{
}

#include "wavelib/wavelib.h"
void tc_compress1(s16 *in, S_PCP_WAVELET &out) //in: 输入, out: 输出缓存
{
	int i=0;
	wave_object obj;
	wt_object wt;
	double *inp,*op;
	inp = new double[S_N];
	op = new double[S_N];
	for (i = 0; i < S_N; ++i)
	{
		inp[i] = in[i];
	}
	int J=3;
	char *name = (char *)"db4";
	obj = wave_init(name);// Initialize the wavelet
	wt = wt_init(obj, "dwt", S_N, J);// Initialize the wavelet transform object
	setDWTExtension(wt, "sym");// Options are "per" and "sym". Symmetric is the default option
	setWTConv(wt, "direct");

	u32 t0 = com_time_getus();
	dwt(wt, inp);// Perform DWT
	u32 t1 = com_time_getus();
	printf("dwt: %d\n",t1-t0);

	for(i=0;i<S_N;i++) //将高频抹去，低频取整
	{
		if(i<ARRAY_SIZE(out.wavelet_buf))
		{
			out.wavelet_buf[i]=round(wt->output[i]);
			wt->output[i] = out.wavelet_buf[i];
		}
		else wt->output[i] = 0;
	}
	//反变换，看误差

	t0 = com_time_getus();
	idwt(wt, inp);
	t1 = com_time_getus();
	printf("idwt: %d\n",t1-t0);

	for(i=0;i<S_N;i++) //与输入的做差
	{
		s16_s_in[i]=round(inp[i]); //主峰时域波形
		//s16_s_in[i]=in[i]-s16_s_in[i]; //与输入做差
		out.dd_list[i] = in[i]-s16_s_in[i]; //与输入做差
	}
	//for(i=1;i<S_N;i++) //误差的差分值
	//{
		//out.dd_list[i] = (out.dd_list[i] - out.dd_list[i-1]);
	//}

	for(i=0;i<S_N;i++) //
	{
		int s16_wave=i<ARRAY_SIZE(out.wavelet_buf)?out.wavelet_buf[i]:0;
		//printf("%d	%d	%d\n", in[i], s16_wave, out.dd_list[i]);
	}
	delete [] inp;
	delete [] op;
}
